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A  class  of  vector-space  bases  is  introduced  for  the  sparse  representation  of  discretizations  of 
integral  operators.  An  operator  with  a  smooth,  non-oscillatory  kernel  possessing  a  finite  number 
of  singularities  in  each  row  or  column  is  represented  in  these  bases  as  a  sparse  matrix,  to  high 
precision.  A  method  is  presented  that  employs  these  bases  for  the  numerical  solution  of  second- 
kind  integral  equations  in  time  bounded  by  <3(nlog2n),  where  n  is  the  number  of  points  in  the 
discretization.  Numerical  results  are  given  which  demonstrate  the  effectiveness  of  the  approach, 
and  several  generalizations  and  applications  of  the  method  are  discussed. 
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1  Introduction 


Integral  equations  are  a  well-known  mathematical  tool  for  formulating  physical 
problems.  As  a  numerical  tool  they  have  several  strengths  (good  conditioning, 
dimensionality  reduction,  and  the  ability  to  treat  arbitrary  regions)  but  have 
had  one  overriding  drawback:  the  high  cost  of  working  with  the  associated  dense 
matrices.  For  a  problem  requiring  an  n-point  discretization,  the  inverse  of  a 
dense  n  X  n-matrix  must  be  applied  to  a  vector.  Even  to  apply  the  matrix 
itself  to  a  vector  requires  order  0(n2)  operations;  application  of  its  inverse  by 
a  direct  (non-iterative)  method  requires  order  0(n3)  operations.  If  an  iterative 
method  is  employed,  the  number  of  iterations  depends  on  the  condition  number 
of  the  problem  and  each  iteration  requires  application  of  the  n  x  n  matrix.  For 
large-scale  problems,  the  resulting  costs  are  often  prohibitive. 

In  recent  years  a  number  of  algorithms  ([5],  [10],  [11],  [14])  have  been  devel¬ 
oped  for  the  fast  application  of  linear  operators  naturally  expressible  as  dense 
matrices,  the  best  known  of  which  are  the  particle  simulation  algorithms  devel¬ 
oped  by  L.  Greengard  and  V.  Rokhlin  [10].  These  schemes  combine  low-order 
polynomial  interpolation  of  the  function  which  defines  the  matrix  elements  with 
a  divide-and-conquer  strategy.  They  achieve  (the  equivalent  of)  order  O(n)  ap¬ 
plication  of  a  dense  n  x  n-matrix  to  a  vector. 

Over  a  somewhat  longer  period,  mathematical  bases  have  been  constructed 
with  certain  multiscale  properties.  Families  of  functions  /?„.{,, 

ha,b{x)  =  \a\~1/2  h^—^j  ,  a, bell,  a  #  0, 

derived  from  a  single  function  h  by  dilation  and  translation,  which  form  a  basis 
for  £2{TZ ),  are  known  as  wavelets  (Grossman  and  Morlet  [12]).  These  families 
have  received  study  by  many  authors,  resulting  in  constructions  with  a  variety  of 
properties.  Meyer  [13]  constructed  orthonormal  wavelets  for  which  h  £  CX('R). 
Daubechies  [7]  constructed  compactly  suppoited  wavelets  with  h  £  C'k( Tv.)  for 
arbitrary  k,  and  [7]  gives  an  overview  and  synthesis  of  the  field. 

A  recent  paper  [6]  establishes  a  connection  between  the  fast  numerical  al¬ 
gorithms  and  the  multiscale  bases.  It  introduces  the  use  of  wavelets  for  the 
application  of  an  integral  operator  to  a  function  in  0(7?  log  n)  operations,  where 
n  is  the  number  of  points  in  the  discretization  of  the  function.  The  dissertation  [4] 
of  one  of  the  present  authors  gives  an  earlier  report  of  the  present  work.  Another 
paper  [2]  constructs  a  class  of  simple  wavelet-like  bases  for  £2[0.  1]  in  which  a 
variety  of  integral  operators  are  sparse.  In  the  present  paper,  rather  than  employ 
a  wavelet  basis  for  £2,  we  construct  a  class  of  bases  that  transform  the  dense 
matrices  resulting  from  the  discretization  of  second-kind  integral  equations  into 
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sparse  matrices.  The  n  X  n-matrices  resulting  from  an  rc-point  discretization  arc 
transformed  to  matrices  with  order  0(n  log  ?i)  non-zero  elements  (to  arbitrary 
finite  precision).  In  these  bases,  the  inverse  matrices  are  also  sparse,  and  are 
obtained  in  order  0(nlog2  n )  operations  by  a  classical  iterative  method  (due  to 
Schulz). 

In  §2  we  present  the  mathematical  construction  of  the  new  bases.  In  §3  we 
briefly  introduce  Nystrom’s  method  for  the  solution  of  integral  equations,  and 
show  how  the  wavelet  bases  result  in  sparse  representation  of  integral  operators 
and  their  inverses.  We  demonstrate  that  the  Schulz  method  of  matrix  inversion 
is  efficient  in  this  context.  In  §4  we  present  the  numerical  algorithms  for  compu¬ 
tation  of  the  wavelet  bases,  transformation  of  an  integral  operator  into  the  bases, 
and  computation  of  its  inverse,  and  we  analyze  the  time  complexity  of  these  algo¬ 
rithms.  A  variety  of  numerical  examples  are  presented  in  §5  to  demonstrate  the 
effectiveness  of  the  approach,  and  generalizations  and  applications  are  discussed 
in  §6. 

2  Wavelet  Bases 

2.1  Properties  of  the  Bases 

Given  a  set  of  n  distinct  points  S  =  {xt,  x2, . . . , xn}  C  71  (the  discretization)  we 
construct  an  orthonormal  basis  for  the  n-dimensional  space  of  functions  defined 
on  S.  For  simplicity,  we  assume  that  n  =  k-2l,  whpre  k  and  1  are  positive  integers, 
and  that  X\  <  x2  <  ■  ■  ■  <  xn.  The  basis  has  two  fundamental  properties: 

1.  All  but  k  basis  vectors  have  k  vanishing  moments,  and 

2.  The  basis  vectors  are  non-zero  on  different  scales. 

Fig.  1  illustrates  a  matrix  of  basis  vectors  for  n  =  128  and  k  =  4.  Each  row  repre¬ 
sents  one  basis  vector,  with  the  dots  depicting  non-zero  elements.  The  first  k  basis 
vectors  are  non-zero  on  xj, . . .  ,x2k,  the  next  k  are  non-zero  on  x2k+i, ....  x.^..  and 
so  forth.  In  all,  one  half  of  the  basis  vectors  are  non-zero  on  2k  points  from  5. 
one  fourth  are  non-zero  on  4 k  points,  one  eighth  are  non-zero  on  8 k  points,  etc. 
Each  of  these  n/2  +  n/4  +  n/S  +  •  •  •  +  =  n  —  k  basis  vectors  has  k  zero  moments, 

be.,  if  6  =  (&i,..., bn)  is  one  of  these  vectors,  then 

n 

]C  b,  ■  -t,3  =  0,  >  =  0,1 . k  -  1. 

i=i 

The  final  k  vectors  result  from  orthogonalization  of  the  moments  (.r1J.  r2J . r„- ) 

bir  7=0  ]  .  I-  —  1 . 


These  properties  of  local  support  and  vanishing  moments  lead  to  efficient 
representation  of  functions  which  are  smooth  except  at  a  finite  set  of  singularities. 
The  projection  of  such  a  function  on  an  element  of  this  basis  will  be  negligible 
unless  the  element  is  non-zero  near  one  of  the  singularities.  As  a  simple  example, 
we  consider  the  function  f(x)  =  log(x)  on  the  interval  [0,1]  with  the  uniform 
discretization  x,  =  i[n.  A  hand  calculation  shows  that  for  any  c  >  0,  /  may  be 
interpolated  on  the  interval  [c,  2c]  by  a  polynomial  of  degree  7  with  error  bounded 
by  4-9,  or  roughly  single  precision  accuracy.  If  we  choose  k  =  8  in  constructing 
the  basis,  /  will  be  represented  to  this  accuracy  by  the  k  basis  vectors  non-zero 
on  ij, . . . ,  X2jt,  the  k  basis  vectors  non-zero  on  Xj, . . . ,  x4fc,  and  so  forth,  down  to 
the  k  basis  vectors  non-zero  on  xi,...,x„,  in  addition  to  the  k  orthogonalized 
moment  vectors.  The  number  of  non-negligible  coefficients  in  the  expansion  of  / 
in  this  basis  grows  logarithmically  in  n,  the  number  of  points  of  the  discretization. 
Although  this  example  is  idealized,  its  behavior  is  representative  of  the  general 
behavior  of  an  analytic  function  near  a  singularity. 


AccuSir/)  T 
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Figure  1:  Hit  matrix  represents  a  wavelet  basis  for  a  discretization  with  128 
points,  for  k  —  4.  Each  row  denotes  one  basis  sector,  with  the  dots  depicting 
non-zero  elements.  All  but  the  final  k  rows  have  k  vanishing  moments. 
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2.2  Construction  of  the  Bases 

The  conditions  of  “local”  support  and  zero  moments  determine  the  basis  vectors 
uniquely  (up  to  sign)  if  vve  require  somewhat  more  moments  to  vanish.  Namely, 
out  of  the  k  vectors  non-zero  onij,...,  £2  k,  we  require  that  one  have  k  vanishing 
moments,  a  second  have  k  +  1,  a  third  have  k  +  2,  and  so  forth,  and  the  k th  have 
2k  —  1  vanishing  moments.  We  place  the  same  condition  on  the  k  basis  vectors 
non-zero  on  x2k+ i,...,X4fc,  and  so  on,  for  each  block  of  k  basis  vectors  among 
the  n  —  k  basis  vectors  with  zero  moments. 

We  construct  the  basis  by  construction  of  a  finite  sequence  of  bases  (shown  in 
Fig.  2),  each  obtained  by  a  number  of  orthogonalizations.  The  first  basis  results 
from  n/(2k)  Gram-Schmidt  orthogonalizations  of  2k  vectors  each.  In  particular, 
the  vectors  (x^, . . . ,  x2k])  for  j  —  0, ...  ,2k  —  1  are  orthogonalized,  the  vectors 
{x2k+\*,  ■  ■  ■ ,  xAk3)  for  j  =  0, . . . ,  2k  —  1  are  orthogonalized,  and  so  forth,  up  to  the 
vectors  (xn^2k+\1  xn])  for  j  =  0, . . . ,  2k  —  1  which  are  orthogonalized. 

Half  of  the  n  vectors  of  the  first  basis  have  at  least  k  zero  moments:  in 
forming  the  second  basis,  these  vectors  are  retained;  the  remaining  n/2  basis 
vectors  are  transformed  by  an  orthogonal  transformation  into  basis  vectors,  each 
of  which  is  non-zero  on  4 k  of  the  points  £i,...,xn,  and  half  of  which  have  at 
least  k  vanishing  moments.  The  orthogonal  transformation  results  from  n/(Ak) 
Gram-Schmidt  orthogonalizations  of  2k  vectors  each.  Similarly,  the  third  basis 
is  obtained  from  the  second  basis  by  an  orthogonal  transformation  that  itself 
results  from  n/(8k)  Gram-Schmidt  orthogonalizations  of  2k  vectors  each.  Before 
we  can  specify  these  orthogonalizations,  we  require  some  additional  notation. 

Suppose  V  is  a  matrix  whose  columns  t’i , . . . ,  v2k  are  linearly  independent .  \\  e 
define  W=Orth(V)  to  be  the  matrix  which  results  from  the  column-by-column 
Gram-Schmidt  orthogonalization  of  V.  Namely,  denoting  the  columns  of  H  by 
u?i, . . . ,  w2k,  we  have 

linear  spanfuq, ...,  u;,}  =  linear  span{iq, ...,  r>,}  .  .  0, 

Wi  10  J  =  dij 

For  a  2k  x  2^-matrix  V  we  let  Vu  and  VL  denote  two  k  x  2A--mat rices.  I  1 
consisting  of  the  upper  k  rows  and  VL  the  lower  k  rows  of  V. 


Now  we  proceed  to  the  definition  of  the  basis  matrices.  Given  the  set  of  points 
Q  =  {'  !■  .,  xnj  C  with  £j  <  •  <  jn,  where  n  —  k  ■  2‘.  we  define  the  2k  x  2k 
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moments  matrices  Mj  ,  for  i 


n /(2k)  by  th^  formula 


Figure  2:  Each  of  the  four  matrices  represents  one  basis,  as  in  Fig.  1.  The 
upper-left  matrix  is  formed  by  orthogonalizing  moment  vectors  on  blocks  of  '2k 
points.  The  upper-right  matrix  is  obtained  from  the  upper-left  matrix  by  pn.mul 
tiplying  by  an  orthogonal  matrix  which  is  the  identity  on  the  upper  half.  Simi 
larly,  the  lower  matrices  are  obtained  by  further  orthogonal  transformations.  Tin 
lo "vr. right  matrix  represents  the  wavelet  basis  for  n  =  6-1,  k  =  4. 


where  5,  =  (i  —  1)2 k.  The  first  basis  matrix  U\  is  given  by  the  formula 


where  U\,tT  =  Orth(A/1,t)  and  nx  =  n/(2k).  The  second  basis  matrix  is  TV-'i , 
with  Ui  defined  by  the  formula 


where  Im  is  the  m  x  m  identity  matrix  and  U'2  is  given  by  the  formula 


where  n2  =  n){Ak ),  U2,,T  =  Orth(  A/2>,-),  and  A U,i  is  given  bv 


M 


2,i  — 


U\,2i-\L  AI\t2i-\  \ 

U\.2iU Ah, 2i  )  ' 


In  general,  the  jih  basis  matrix,  for  j  =  2, . . . ,  \og2(n/k).  is  with! 

defined  by  the  formula 

h  —  n/2J~i 

f; 
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where  U'  is  given  by  the  formula 


(  u* 


Uj.1 


u 


U*L 


TJ.  v 

UJ,2 


Ui.n,' 


\  Uj.n,V  ) 

where  n}  =  nj{23k),  Uh,  is  given  by 

Uh,T  =  Orth(A/,-,j), 

and  Mhl  is  given  by 


Mh>  -  TT  u  ,f 


(2) 


(3) 


The  final  basis  matrix  U  =  Ur  ■  ■  U\,  where  /  =  log 2(?i/t),  represents  the  wavelet 
basis  of  parameter  k  on  X\ , . . . ,  xn. 


Remark  2.1  The  definitions  given  for  the  basis  matrices  are  mathematical  def¬ 
initions  only;  in  a  numerical  procedure,  considerable  roundoff  error  would  be 
introduced  by  the  orthogonalizations  defined  above.  In  the  actual  implementa¬ 
tion,  the  matrices  MJtX  are  shifted  and  scaled,  resulting  in  a  numerically  stable 
procedure  that  is  equivalent  to  the  above  definitions  (in  exact  arithmetic).  Details 
of  this  procedure  are  provided  in  §4. 


It  is  apparent  that  the  application  of  the  matrix  V  to  an  arbitrary  vector 
of  length  n  may  be  accomplished  in  order  0(n )  operations  by  the  application 
of  U i,. . .  ,17/  in  turn.  Similarly,  U~l  =  UT  may  be  applied  to  a  vector  in  order 
O(n)  operations.  Certain  dense  matrices,  in  particular  those  arising  from  integral 
operators,  are  sparse  in  the  basis  of  U  and  their  similarity  transformations  can 
be  computed  in  0(n  log  n)  operations.  These  techniques  are  developed  in  the 
following  sections. 

We  conclude  this  section  with  an  illustration  of  the  vectors  of  one  basis  from 
this  class,  in  Fig.  3. 
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Figure  3:  Basis  vectors  on  four  scales  are  shown  for  the  basis  where  n  =  128. 
points  Xi,...,xn  are  equispaced,  and  k  =  8.  The  first  column  of  vectors  consists 
of  rows  1-8  of  U,  the  second  column  consists  of  rows  65-72,  etc.  .Vo/e  that  halj 
of  the  vectors  are  odd  and  half  are  even  functions,  and  that  the  odd  ones  an 
generally  discontinuous  at  their  center. 
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3  Second-Kind  Integral  Equations 
3.1  Nystrom  Method 

A  linear  Fredholm  integral  equation  of  the  second  kind  is  an  expression  of  the 
form 

/(i)  -  f  K(x,t)  f(t)  dt  =  g(x),  (4) 

J  a 

where  the  kernel  I\  is  in  C2[a,  b j2  and  the  unknown  /  and  right-hand-side  <j  are  in 
£2[a,6],  We  use  the  symbol  K.  to  denote  the  integral  operator  of  Eq.  (4).  which 
is  given  by  the  formula 

(A?/)(x)  =  /  J\(x,t)  fit)  dt, 

J  a 

for  all  /  £  £2[a,6]  and  x  £  [a,b],  Then  Eq.  (4)  written  in  operator  form  is 

{1-K)f  =  cj.  (5) 

The  Nystrom,  or  quadrature,  method  for  the  numerical  solution  of  integral  equa¬ 
tions  approximates  the  integral  operator  K  by  the  finite-dimensional  operator  ft. 

characterized  by  points  Xi,i2, . . .  ,x„  €  [a,b]  and  weights  U'j.ug . w„  £  R. 

and  given  by  the  formula 

(Rf)(x)  =  jr  Wj  I\{x,Xj) 

J=1 

for  all  /  £  £2[a.6J  and  x  £  [a,6J.  Substitution  of  R  for  JC  in  Eq.  (5).  combined 

with  the  requirement  that  the  resulting  equation  hold  for  x  =  Xi,x2 . xn. 

yields  the  following  system  of  n  equ;  .ions  in  the  n  unknowns  /j,  /2, . .  . ,  fn: 

n 

/.  -  lL’:  A'(x’t. x-j )  fj  -  g(x, J.  i  =  1 . V.  ((»! 

r=i 

The  approximation  (/i,. .  .  ,  fn)  to  the  solution  /  of  Eq.  (4)  may  be  extended 
to  all  x  £  [a,  6]  by  the  natural  formula 

n 

//?(■*•)  =  <7(d  +  Ylu'j  R(x.Xj)  /}.  ("I 

1=1 

which  satisfies  fn{x,)  =  f,  for  i  —  \  , ,r>.  How  large  is  the  error  c. /?  =  /  —  j),-  i<! 
the  approximate  solution?  We  follow  the  derivation  by  Delves  and  Mohamed  in 
[S].  Rewriting  Eqs.  (6)  in  operator  form,  we  have 

(/  -  R)fR  =  g.  '  v  1 
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and  combining  Eqs.  (5)  and  (S)  yields 


(/  -  K,)cr  =  (AC  —  R)/r. 

Provided  that  (I  —  AC)-1  exists,  we  obtain  the  error  bound 

l|e«||  <  ||(/ - AT)-1!!  -  ||(AC  -  (U) 

The  error  depends,  therefore,  cn  the  conditioning  of  the  original  integral  equation, 
as  is  apparent  from  the  term  ||(7  —  AC)-1||,  and  on  the  fidelity  of  the  quadrature 
R  to  the  integral  operator  AC.  It  is  not  necessary  that  l|AC  —  i?||  be  small,  rather 
merely  that  R  approximate  AC  well  near  the  solution  /.  Quadrature  rules  that 
have  this  property,  but  which  j.re  defined  only  on  the  points  Xj,  .  .  . ,  xn,  are  de¬ 
veloped  in  [3].  In  these  rules  the  quadrature  weights  wX]  depend  on  the  point 
x,,  for  i  =  l,...,n,  and  the  quadrature  rules  converge  rapidly  for  kernels  with 
diagonal  singularities.  These  rules  are  used  below  in  the  numerical  examples  of 

§5- 


3.2  Sparse  Representation  of  Integral  Operators 

We  concern  ourselves  here  with  kernels  K  =  I\(x,t)  which  are  analytic  except 
at  x  =  t,  where  they  possess  an  integrable  singularity.  We  initially  discretize 
the  integral  operator  AC  using  a  simple  equispaced  quadrature.  Given  n  >  2.  we 
define  points  Xi,....x„  to  be  equispaced  on  the  interval 

x,  =  a  +  (i  -  1)(6  -  a)/(n  -  1),  (10) 

and  define  the  elements  Tt]  of  the  n  x  n-matrix  T  by  the  formula 

T  =  S  '*i  nn 

u  \  0  i  =  j. 

Note  that  the  matrix  T  —  T(n)  corresponds  to  a  primitive,  trapezoid-like  quadra¬ 
ture  discretization  of  the  integral  operator  AC.  The  matrix  T  possesses  the  same 
smoothness  properties  as  the  kernel  A'(x,<).  Transformation  of  T  by  the  bases 
developed  in  §2  produces  a  matrix  that  is  sparse,  to  high  precision.  The  number 
of  elements  is  effectively  bounded  by  order  0{n  log  ) . 

When  the  matrix  representing  the  quadrature  corrections  developed  in  [31 
is  added  to  T,  producing  high-order  convergence  to  the  integral  operator,  this 
complexity  bound  remains  valid. 

The  matrix  T,  transformed  by  the  orthogonal  n  x  n-matrix  U ,  can  be  de¬ 
composed  into  the  sum  of  a  sparse  matrix  and  a  matrix  with  small  norm.  Given 
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t  >  9,  th^re  exists  cf  >  0,  independent  of  n,  such  that  the  transformed  matrix- 
can  be  written  in  the  form 

UTUt  =  V  +  E, 

where  the  number  of  elements  in  V  ~  V(n)  is  bounded  by  ct  n  log  n  and  E  =  E(n) 
is  small:  Hit’ll  <  ej|J’||.  We  do  not  prove  this  assertion  here;  the  proof  parallels 
the  proofs  of  similar  statements  in  [2],  but  is  somewhat  more  tedious. 


3.3  Solution  via  Schulz  Method 


The  sparse  matrix  representing  the  integral  operator  also  has  a  sparse  inverse, 
which  can  be  computed  rapidly. 

Schulz’s  method  [15]  is  an  iterative,  quadratically  convergent  algorithm  for 
computing  the  inverse  of  a  matrix.  Its  performance  is  characterized  as  follows. 


Lemma  3.1  Suppose  that  A  is  an  invertible  matrix,  A’o  is  the  matrix  given  by 
X0  =  Ah/\\AhA\\  ,  and  form  =  0,1,2,...  the  matrix  Xm+l  is  defined  by  the 
recursion 


X 


m  +  1 


=  2Xm  -  XmAX, 


Then  Xm+\  satisfies  the  formula 


I  -  Xm+\A  -  (I  -  XmA)2. 


(12) 


Furthermore,  Xm  — >  -4-1  as  m  —>  oc  and  for  any  e  >  0  we  have 


|| /  -  .Vm/4||  <  c  provided  m  >  2  log2  k(A)  +  log2  log(l/e),  (13) 


wheie  k(.4)  =  ||.4||  •  1 1 .4 ~ 1  j ;  is  the  condition  number  of  A  and  the  norm  is  yirin 
by  ||/1||  =  (largest  eigenvalue  of  A11  Afil2 . 


Proof  Eq.  (12)  is  obtained  directly  from  the  definition  of  .Vm+i-  Bound  (13) 
is  equally  straightforward.  Noting  that  Alf  A  is  symmetric  positive-definite  and 
letting  A0  denote  the  smal'est  and  \l  the  largest  eigenvalue  of  ,4;/.4  we  have 


/  -  A  o-4 1 1 


A11 
II -4". 4 


i  —  -W-i] 

1  -  s-(.4)-\ 


(13) 


From  Eq.  (12)  we  obtain  /  —  Am,4  —  (I  —  A0.4)2m.  which  in  combination  with 
Eq.  (14)  and  simple  manipulation  yields  bound  (13).  □ 

The  Schulz  method  is  a  notably  simple  scheme  for  matrix  inversion  and  it' 
convergence  is  extremely  rapid.  It  is  rarely  used,  however,  because  it  involves 


11 


matrix-matrix  multiplications  on  each  iteration;  for  most  problem  formulations, 
this  process  requires  order  0(n3)  operations  for  an  n  x  n  matrix.  We  observe, 
however,  that  a  sparse  matrix,  possessing  a  sparse  inverse,  whose  iterates  Xn 
are  also  sparse,  may  be  rapidly  inverted  using  the  Schulz  method.  We  have 
seen  above  that  a  discretized  integral  operator  I  —  T,  similarity-transformed  to 
the  representation  A  =  I  —  UTU T,  has  only  order  O(nlogn)  elements  (to  finite 
precision).  In  addition,  AT A  and  ( ATA)m  are  similarly  sparse.  This  property 
enables  us  to  employ  the  Schulz  algorithm  to  compute  A~x  in  order  0{n  log"  n) 
operations. 

3.4  Oscillatory  Coefficients 

We  now  consider  a  somewhat  more  general  class  of  integral  equations,  in  which 
the  integral  operator  is  given  by  the  formula 

(D£f)(x)  =  p(x)  [b 

J  a 

where  the  kernel  K  is  assumed  to  be  smooth,  but  the  coefficient  function  p 
can  be  oscillatory.  In  particular,  we  only  restrict  p  to  be  positive.  In  terms  of 
generality,  these  problems  lie  between  the  problems  with  smooth  kernels  (and 
constant  coefficient)  and  those  with  arbitrary  oscillatory  kernels. 

Writing  the  corresponding  integral  equation  in  operator  form,  we  obtain  the 
equation 

(/  —  DK)f  =  g.  (15) 

Although  D  is  a  diagonal  operator,  and  K.  is  smooth,  it  is  clear  that  the  discretiza¬ 
tion  of  the  operator  DK  will  not  be  a  sparse  matrix  in  wavelet  coordinates.  In 
this  framework,  it  would  appear  that  the  construction  of  this  paper  is  inapplica¬ 
ble.  If  we  instead  consider  the  operator  DX^2KDX^2 ,  in  which  oscillations  in  the 
rows  match  those  in  the  columns,  it  becomes  clear  that  the  construction  of  §2 
can  be  revised.  Rather  than  constructing  basis  functions  orthogonal  to  low-order 
polynomials  x: ,  we  can  construct  them  to  be  orthogonal  to  p{x)1^2  xJ.  The  sole 
revision  in  our  definition  of  basis  matrices  Ui, . . .  ,Ui  is  to  replace  the  definition 
( 1 )  of  the  moments  matrices  A/I  (  for  i  =  1, . . . ,  nj(‘2k),  by  the  new  definition 


t  Ps,  + 1 

0 

0  \ 

(  1  JT+1  •' 

T  2fc-1 

xs,  +  l 

\ 

Mu  = 

0 

Ps,+7 

0 

1  W,  +  2  •  ' 

x  2k—  1 

xs,+2 

\  0 

0  Ps,  +  2k  / 

\  1  Xs,  +  ?k  ■■ 

r  ,2k- 1 

xs,+7k 

/ 

where  s,  =  (i  —  1  )2A-  and  p}  =  p(j;)1/2. 
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Now  the  integral  equation  (15)  can  be  transformed  to  the  equation 
(/  -  D1/2fCD1/2)(D~1/2f)  =  (D~1/2g), 

which  is  discretized  to  a  system  that  is  sparse  in  the  revised  wavelet  coordinates. 
The  inverse  matrix  is  also  sparse. 

4  Numerical  Algorithms 

In  §2  we  defined  a  class  of  bases  for  functions  defined  on  {ij, . . . ,  xn}  and  in  §3  we 
showed  that,  to  finite  precision,  second-kind  integral  operators  and  their  inverses 
are  asymptotically  sparse  in  these  bases.  In  this  section  we  present  procedures 
for  computation  of  the  bases,  discretized  integral  operators  in  these  bases,  and 
the  inverses  of  these  operators.  In  §5  we  give  some  numerical  examples  based  on 
our  implementations  of  these  procedures. 

The  computation  of  the  wavelet  bases  is  discussed  next,  followed  by  a  dis¬ 
cussion  of  the  transformation  of  the  integral  operators  to  the  wavelet  bases.  We 
defer  discussion  of  the  computation  of  the  inverses,  sketched  above,  to  §4.3, 
which  contains  detailed  descriptions  of  all  of  the  algorithms.  Finally,  §4.4  gives 
the  complexity  analysis  for  the  algorithms. 


4.1  Computation  of  Wavelet  Bases 


It  was  mentioned  in  §2  that  the  mathematical  definition  of  if  used 

directly,  would  result  in  a  numerical  procedure  that  would  create  large  roundoff 
errors.  A  correct  procedure  is  obtained  by  shifting  and  scaling  the  matrices  3/,., 
defined  there. 

For  a  pair  of  numbers  (p,  cr)  6  TZ  x  (7£\{0} )  we  define  a  2k  x  2&-matrix  S(/i.'t) 
whose  (i.  j)th  element  is  the  binomial  term 


(j-i) 

\i-  l)  <j;-' 


(10) 


for  i  <  j,  and  5(p,<7),iJ  =  0  otherwise.  The  matrix  S(//,< 7)  is  upper-triangular 
and  non-singular,  and  its  inverse  is  given  by  the  formula 

S(p.cr)"1  =  5(-p/cr,  1/a).  (IT) 


Furthermore,  the  product  formula 


S(n  1 .  <7\  Q-l )  —  S(fl  1  +  Hl<J I  .  <71  <7-1 ) 


US' 


is  easilv  verified. 
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We  define  A/j,-  for  j  =  1, . . . ,  /  and  i  =  1, . . . ,  n/(2Jk)  by  the  formula 

A/',  =  A^.S^,,,^,),  (19) 

where  fiJti  =  (xl+li_l)k2,  +  xtk2,)/2,  ah,  =  (xik2,  -  xl+[i_l)kv)/2  and  the  matrix 
A/ji  is  defined  by  Eq.  (1)  and  Eq.  (3)  in  §2.  The  matrix  Uh%  is  given  by  the 
formula 

U:/  =  Orth  (A/',),  (20) 

which  is  equivalent  to  the  definition  given  by  Eq.  (2).  This  equivalence  immedi¬ 
ately  follows  from  the  fact  that  S(n,cr)  is  upper-triangular  and  non-singular. 
The  matrices  A/j  ,  for  i  =  1, . . .  ,n/(2£)  are  actually  computed  by  the  formula 


Mu  = 


/  J  +  1  -Ml,. 


/ 1 3,  -f  2  —  Ml  ,i  j  2k~1 


X3,+2k-Vl,> 


2k—  1 


/ 


where  s,  =  (i  —  1)2E.  Likewise,  the  matrices  A/j,  for  j  =  2, 
1, . . . ,  n/(2Jk)  are  computed  by  the  formula 

Vi 


M'  _  [  M']-\,2x-\S},t 


where  Sj,  and  Sj,  are  defined  by  the  formulae 

Syj  S(/ij_ l,2t—  1,  &]  —  1 ,2*  —  1  )  S(/iy,,  ffy,  ) 

Sjtt  —  S  &j-l,2i)  (Jjj) . 


(21) 

,  /  and  i  = 

(22) 


(23) 

(24) 


Application  of  the  inverse  and  product  rules  given  in  Eq.  (17)  and  Eq.  (IS)  to 
Eq  (23)  and  Eq.  (24)  yields  formulae  by  which  Sj,  and  Sj,  can  be  computed: 

Sj,  S( (/iyi  pj  — 1,2«—  1  )/Sj  —  1 ,2i  —  1  i  0j,i [Vj  —  1 ,2<—  1 )  (-'^) 

Sj,  =  S((/ij,,  —  ^>  —  1,2*-.  < Tj,i/(TJ-I,2i )■  (2b) 

The  matrices  A/j,  given  by  Eq.  (21)  and  Eq.  (22)  are  easily  seen  to  be  math¬ 
ematically  equivalent  to  those  defined  by  Eq.  (19);  nonetheless,  computation  of 
A/j,  using  Eqs.  (21)  and  (22)  avoids  the  large  roundoff  errors  which  would  oth¬ 
erwise  result. 


4.2  Transformation  to  Wavelet  Bases 

We  assume  that  for  equispaced  points  (defined  in  Eq.  (10))  and  some 

k  the  orthogonal  matrices  defined  in  §2  have  been  computed  (/  = 

log2(n/k)).  We  now  present  a  procedure  for  computation  of  UTUt .  where  (  = 
L'l  ■  ■  ■  i'i  and  T  is  the  discretized  integral  operator  defined  in  Eq.  (11). 
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4.2.1  Simple  Example 

We  begin  with  a  simplified  example  in  which  T  is  replaced  by  an  n  x  n-matrix 
V  of  rank  k  whose  elements  are  defined  by  the  equation 

k  k 

Vii  —  y  ]  y  (  A raX{  Xj  ,  i,  j  =  1 , . . . ,  n. 

r=l  3=1 

Each  row  and  each  column  of  V  contain  elements  which  are  the  values  of  a 
polynomial  of  degree  k  —  1.  The  matrix  V  can  be  written  as  V  =  PTAP ,  where 
the  elements  of  the  k  x  n-matrix  P  are  defined  by  P,2  =  i/-1  and  A  is  the 
k  x  fc-matrix  with  elements  AtJ.  Recalling  that  the  last  k  rows  of  the  basis 

matrix  U  consist  of  an  orthogonalization  of  the  moment  vectors  (xp . x,p)  for 

j  =  0, —  1,  we  can  rewrite  V  as  V  =  (P')T  A'P'.  Here  the  k  x  n-matrix  P' 
consists  of  the  last  k  rows  of  U  and  A'  is  a  new  k  x  fc-matrix  with  elements  A'-  . 

By  the  orthogonality  of  £/,  it  is  clear  the  n  xn-matrix  UVUT  —  U(P')T A1  P'UT 
consists  entirely  of  zero  elements  except  the  k  x  fc-submatrix  in  the  lower-right 
corner,  which  is  the  matrix  A'.  Given  a  function  to  compute  elements  of  the 
n  x  n-matrix  V ,  the  matrix  A'  can  be  computed  in  time  independent  of  n  by 
using  a  k  x  k  extract  of  values  from  V.  We  form  the  matrix  V'  with  elements  \\: 
defined  by  the  formula 


—  Kn/l-jn/b  PJ  —  l,...,fc.  (-1) 

Then  V'  =  (P")T A' P" ,  where  P"  is  the  k  x  k  extract  of  P'  with  elements  given 
by  P”  =  P- jn/k.  Thus  we  obtain 

A '  ={{P")T)-1V'{P!,)-X  (2S) 

from  P"  and  V'  readily  in  0(k3)  operations,  and  we  have  obtained  UVUT . 

4.2.2  General  Case 

The  integral  operator  matrix  T  is,  of  course,  not  of  low  rank,  but  it  can  be  divided 
into  submatrices,  each  approximately  of  rank  k  (see  Fig.  4).  The  submatrices 
near  the  main  diagonal  are  of  size  k  x  k,  those  next  removed  are  2 k  x  2k  and 
so  forth  up  to  the  largest  submatrices,  of  size  n/ 4  x  n/ 4.  The  total  number  of 
submatrices  is  proportional  to  n/k.  Given  an  error  tolerance  e  >  0.  k  may  be 
chosen  (independently  of  n)  so  that  each  submatrix  of  T.  say  T\  may  be  written 
as  a  sum.  T'  =  V'  +  E' ,  where  the  elements  of  V71  are  given  by  a  polynomial  of 
degree  k  —  1  and  ||£‘||  <  e||P||. 

The  simplified  example,  in  which  the  matrix  to  be  transformed  is  of  rank 
k,  is  now  applicable.  Each  submatrix  of  T  is  treated  as  a  matrix  of  rank  k 
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Figure  4:  The  matrix  represents  a  discretized  integral  operator  with  a  kernel  that 
is  singular  along  the  diagonal.  The  matrix  is  divided  into  submatrices  of  rank  k 
(to  high  precision)  and  transformed  to  a  sparse  matrix  with  O(nlogn)  elements. 
Here  n/k  =  32. 

and  is  transformed  to  wavelet  coordinates  (for  its  own  scale)  in  order  0(k 3) 
operations.  To  make  this  precise,  we  write  T  =  T0  +  •  •  •  +  7}_2  where  T,  consists 
of  the  submatrices  of  size  2'k  x  2lk.  For  each  the  submatrices  of  T,  may  be 
interpolated  by  rank  k  submatrices,  as  indicated  by  the  extract  of  Eq.  (27),  to 
obtain  matrices  V,.  Thus  T,  =  V, ■  +  E, ,  where  j|£',||  is  small.  In  the  simplified 
example  above,  we  have  shown  that  the  transformed  matrices 

W0  =  Vo 

Wl  =  UXVXUXT 

W2  =  U2UxV2UxtU/  (29) 

W/-2  =  Ul-2---UlVi-2U1T---Ul.2T 


can  be  computed  by  many  applications  of  Eq.  (2S),  all  in  order  0(nk 2)  operations. 
This  estimate  follows  from  the  fact  that  there  are  0{n/k)  submatrices,  each 
of  which  is  transformed  in  0[k3)  operations.  Now  we  define  n  x  n-mat rices 


Ro,  ■  ■  ■  ,Ri  recursively: 


(W0  i  =  0 

{  UiRi-iUf  +  Wi  i>  1 


(30) 


(here  j  =  Wi  =  0).  Then  R[  contains  the  final  result,  Ri  =  U(T  —  E)VT , 
where  E  =  E0  +  •  •  •  +  Ei- 2. 

The  matrix-matrix  products  in  the  definition  of  Ro, . . .  ,Ri  can  be  computed 
directly,  since  the  factors  and  the  products  contain  no  more  than  0(n  logn)  ele¬ 
ments.  A  simple  implementation  with  standard  sparse  matrix  structures  results 
in  a  total  operation  count  of  order  0(n  log2  n),  but  an  implementation  using 
somewhat  more  elaborate  data  structures,  in  which  repetitive  handling  of  data 
is  avoided,  requires  only  order  O(nlogn)  operations. 

Computation  using  the  result  Ri  is  made  more  efficient  by  removing  the  ele¬ 
ments  of  Ri  which  can  be  neglected,  within  the  precision  with  which  Ri  approx¬ 
imates  UTUT.  For  a  given  precision  e,  we  discard  a  matrix  E1  by  eliminatin 
elements  from  Rt  below  a  threshold  r.  The  threshold  depends  on  the  choice  ol 
norm;  in  our  implementation,  we  use  the  row-sum  norm 


||  A||  =  max^lA.^I, 

‘  j=i 

for  an  n  x  rc-matrix  A.  The  element  threshold 

r  =  ~\\T\\  (31) 

n 

clearly  results  in  a  discarded  matrix  E'  with  ||£’,||  <  e||T||. 


4.3  Detailed  Descriptions  of  Algorithms 


Procedure  to  compute 

Comment  [Input  to  this  procedure  consists  of  the  number  of  points  ??. 
the  number  of  zero  moments  k,  and  the  points  Xi,...  ,xn.  Output  is  the 
matrices  Uhl  for  j  =  1, . . . ,  /  and  i  =  1, . . . ,  n/(2Jk),  which  make  up  the 
matrices  U\,. . .  ,Ui  (note  /  =  log2(n/^'))-] 

Step  1. 

Compute  the  shifted  and  scaled  moments  matrices  M[ ,  for 
i  =  1, . . .  ,n/(2k)  according  to-  Eq.  (21). 


'cO  ' 


Step  2. 

Compete  £/],,  from  M[  t  by  Eq.  (20)  using  Gram-Schmidt 
orthogonalization  for  i  —  1, . . . ,  n/(2k). 

Step  3. 

Comment  [Compute  A/)  1  and  Uhl  for  j  =  2, . . . ,  /  and 
*'  =  l,...,n/(2>fc).] 
do  j  =  2 

do  i  =  1, . . .  ,n/(23k) 

Compute  and  {/,_ 2>. 

Compute  5],  by  Eq.  (25)  and  S£t-  by  Eq.  (26); 

multiply  to  obtain  A/j,  by  Eq.  (22). 
Orthogonalize  M'}1  to  obtain  Uhi  by  Eq.  (20). 

enddo 

enddo 


Procedure  to  compute  UTUT 

Comment  [Input  to  this  procedure  consists  of  n,  k,  the  matrices  Uh, 
computed  above,  a  function  to  compute  elements  of  T,  and  the  chosen 
precision  e.  Output  is  a  matrix  Ri  such  that  ||/?(  —  UTUT\\  <  e|[ Z’H .] 

Step  4. 

Compute  the  k  x  k  extracts,  indicated  by  Eq.  (27),  of  the 
submatrices  of  T  shown  in  Fig.  4. 

Step  5. 

Extract  the  matrices  P"  (Eq.  (28))  from  C/2C’i,  . . . ,  Ur  ■  -U\ 
and  compute  VE0, ....  H//_2  according  to  Eqs.  (29). 

Step  6. 

Compute  Rq,...,R[  by  Eq.  (30),  discarding  elements  below 
a  threshold  r  determined  by  the  precision  e  (Eq.  (31)). 


Procedure  to  compute  UT  1UT 

Comment  [Input  to  this  procedure  consists  of  n,  the  matrix  /?;  which 
approximates  UTUT  ,  and  the  precision  e.  Output  is  a  matrix  A’m 
that  approximates  UT~lUT .} 
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Step  7. 

Compute  the  matrix  Xo  =  RtT  Ri/\\RiT  Ri\\  by  direct  matrix 
multiplication,  discarding  elements  below  a  threshold  r 
determined  by  the  precision  e  (Eq.  (31)). 

Step  8. 

Comment  [Obtain  the  inverse  by  Schulz  iteration.] 
do  m  =  0, 1, . . .  while  |[7  -  XmRi\\  -  t 

Compute  Xm+\  =  2Xm  —  XmRiXm,  discarding  elements 
below  threshold, 
enddo 

4.4  Complexity  Analysis 

In  the  following  table,  we  provide  the  operation  count  for  each  step  of  the  com¬ 
putation  of  UT~lUT . 

Step  Complexity  Explanation _ 

1.  0(nk)  There  are  n/(2k)  2k  x  2fc-matrices;  each  element 

of  the  matrices  is  computed  in  constant  time. 

2.  0(nk2)  For  each  of  the  n/(2k)  matrices,  perform  a  Gram- 

Schmidt  orthogonalization  requiring  order  0(k 3) 
operations. 

3.  0(nk2)  For  each  of  n/(4fc)  -f  n/(8k)  +  •  •  •  +  1  =  n/(2k)  —  1 

matrices,  compute  four  products  of  a  k  x  2k- 
matrix  with  a  2k  x  2fc-matrix,  construct  two 
2k  x  2A:-matrices,  and  orthogonalize  one  2k  x  2k- 
rnatrix. 

4.  O(nk)  There  are  6(1 +3 +7 -| - f  (n/(2fc)  — 1))+ 3(n/A’)  — 

2,  or  order  0(n/k),  submatrices  of  T  and  for  each 
matrix  we  compute  k2  elements. 

5.  0{nk2)  There  are  n/(2A:)  +  n/(4A')  4 - \- 1  =  n/ k  —  1  ma¬ 

trices  P" ,  each  the  product  of  two  k  x  ^-matrices. 

These  are  each  inverted  and  multiplied  with  the 
0{n/k)  matrices  of  the  previous  step. 
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Step  Complexity  Explanation 

6.  0(n  log  n)  The  diagonally-banded  matrix  Wo,  which  contains 

0{n)  elements,  grows  to  O(nlogn)  elements  by 
the  computation  of  UW0UT ,  as  can  be  seen  by 
simply  examining  pictures  of  W0  and  U.  The  non¬ 
zero  elements  of  the  transformed  \\\ , . . . ,  Wj_2  are 
a  subset  of  those  of  Wo. 

7.  0(nlog2n)  Multiplication  of  two  matrices,  each  with  order 

O(nlogn)  elements,  to  obtain  a  product  with  or¬ 
der  O(nlogn)  elements. 

8.  0(nlog2n)  Two  multiplications  like  that  of  Step  7  are  made 

per  iteration;  the  number  of  iterations  is  indepen¬ 
dent  of  n  and  given  by  bound  (13). 

Total  0(n\og2n) 

5  Numerical  Examples 

In  this  section  we  present  operators  from  several  integral  equations,  the  discretiza¬ 
tion  and  transformation  of  the  operators  to  our  wavelet  bases,  and  the  inversion 
of  the  operators  via  Schulz  method. 

5.1  Uncorrected  Quadratures 

We  first  examine  simple  quadratures  with  equal  weights,  except  weight  zero  at 
the  singularity,  as  represented  by  matrix  T  =  T(n)  defined  by  Eq.  (11).  We 
transform  the  matrix  I  —  T  to  wavelet  coordinates  as  described  in  §4.2,  then 
compute  (I  —  T)~l. 

These  discretizations  are  not  particularly  useful  for  the  solution  of  the  inte¬ 
gral  equations,  due  to  their  slow  convergence  to  the  integral  operators.  They 
nonetheless  make  good  illustrative  examples,  for  they  retain  the  smoothness  of 
the  operator  kernels  and  produce  correspondingly  sparse  matrices.  In  the  next 
subsection,  we  examine  the  results  of  using  high-order  quadratures. 

For  various  sizes  n  of  discretization,  we  tabulate  the  average  number  of  ele¬ 
ments  per  row  in  the  transformed  matrix  U{I  —  T)UT  and  the  computation  time 
to  obtain  the  matrix.  In  addition,  we  display  the  average  number  of  elements 
per  row  of  its  inverse,  and  the  time  to  compute  the  inverse.  Finally,  we  show  the 
error  introduced  by  these  computations.  The  error  is  determined  by  the  appli¬ 
cation  of  the  forward  and  inverse  transformations  to  a  random  vector.  Choose  a 


vector  v  of  length  n  with  uniformly  distributed  pseudo-random  elements;  com¬ 
pute  ( I  —  T)v  directly,  by  a  standard  procedure  requiring  order  0(n 2)  operations; 
transform  to  wavelet  coordinates,  obtaining  U{I—T)v\  apply  the  computed  value 
of  if  (I  —  T)~xUt  to  the  vector  U(I  —  T)v,  transform  to  original  coordinates  by 
application  of  UT\  compare  the  result  v'  to  v.  The  measure  of  error  is  the  relative 
C2  error,  defined  by  the  formula 


V  !  1 


The  programs  to  transform  and  invert,  as  well  as  those  to  determine  the  error, 
were  implemented  in  FORTRAN.  All  computations  were  performed  in  double¬ 
precision  arithmetic  on  a  Sun  Sparcstation  1. 

The  first  set  of  examples  is  for  the  kernel  K(x,  t)  =  log  |x  —  t\,  for  a  wavelet 
basis  of  order  k  =  4  and  various  choices  of  precision  t.  The  matrix  sparsities, 
execution  times,  and  errors  appear  in  Table  1.  Although  the  sparse  matrices  are 
not  banded,  we  loosely  refer  to  the  average  number  of  matrix  elements  per  row 
as  the  matrix  bandwidth.  We  make  the  following  observations: 


1.  The  band  widths  N\,N?  of  the  operator  and  its  inverse  decrease  with  increas¬ 
ing  matrix  size.  In  other  words,  in  the  range  of  matrix  sizes  tabulated,  the 
number  of  matrix  elements  grows  sublintarly  in  the  matrix  dimension  n. 


2.  The  operator  matrix  in  wavelet  coordinates  is  computed  in  time  that  grows 
nearly  linearly  in  n. 


3.  The  inverse  matrix  is  computed  in  time  which  grows  sublinearly  in  n.  This 
is  due  to  the  fact  that  the  cost  of  multiplying  the  sparse  matrices  is  roughly 
order  0(nN2),  for  size  n  and  bandwidth  N.  One  result  is  that  the  cost 
sometimes  drops  as  n  increases. 


4.  The  accuracy  is  within  the  precision  specified.  In  fact,  due  to  the  conser¬ 
vative  element  thresholding  (Eq.  31),  the  actual  error  is  considerably  less 
than  e. 


5.  The  cost  increases  with  increasing  precision  e,  due  to  the  increasing  band- 
widths  generated.  The  the  bandwidths  increase  approximately  as  log(l/e). 

6.  For  k  =  4,  our  fast  transformation  algorithm  does  not  maintain  the  specified 
precision  of  t  =  10-4.  This  anticipated  result  follows  from  the  error  estimate 
for  polynomial  interpolation  of  logarithm  on  intervals  separated  from  the 
origin.  An  unanticipated  attendant  result  is  that  the  bandwidth  increases 
as  the  quality  of  approximation  deteriorates  (compare  to  k  =  S.  below).  As 
a  result,  we  did  not  complete  examples  for  ji  =  409G.  8192. 
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Table  1:  The  integral  operator  1C  defined  by  the  formula  (fCf)(x)  =  f(x)  — 
/0!  log  |x  —  <|  /(<)  dt  is  discretized,  transformed  to  the  wavelet  coordinates  with 
k  =  4,  and  inverted.  For  various  precisions  e  and  various  sizes  of  discretization, 
we  tabulate  the  average  number  of  elements/row  N\  of  the  matrix  in  wavelet  co¬ 
ordinates  and  the  time  in  seconds  t j  to  compute  it,  corresponding  statistics  :V2 
and  <2  for  the  inverse,  and  the  error  (see  text). 


Transform.  Inversion  C? 

e  n  N\  1 1'  ftr2  <2  Error 

1(H  64  72  2  O  2  0.503E-02 

128  6.9  3  6.5  4  0.257E— 02 

256  3.3  7  4.4  4  0.250E-02 

512  2.8  13  3.1  6  0.236E-02 

1024  1.0  26  2.1  6  0.227E-02 

2048  1.4  40  1.4  6  0.221E-02 

4006  1.2  07  1.2  8  0.221E-02 


10~4 


8102 

1.1 

105 

1.1 

12 

0.217E-02 

64 

17.6 

2 

10.5 

14 

0.350E-03 

128 

18.1 

5 

20.0 

36 

0.270E-03 

256 

18.0 

11 

•  20.0 

83 

0.331E-03 

512 

14.5 

21 

15.7 

123 

0.257E-03 

1024 

13.3 

41 

15.5 

262 

0.340E— 03 

2048 

8.5 

73 

0.8 

287 

0.233E-03 

4006 

5.8 

131 

6.5 

304 

0.222E-03 

8102 

3.7 

242 

4.4 

312 

0.221E-03 

64 

28.4 

3 

30.3 

36 

0.104E-03 

128 

32.1 

6 

34.3 

111 

0.140E-03 

256 

34.5 

15 

37.5 

302 

0.161E-03 

512 

33.1 

31 

35.8 

618 

0.177E-03 

1024 

30.2 

63 

33.6 

1280 

0.180E-03 

2048 

25.0 

121 

27.6 

2040 

0.102E-03 

7.  The  inversion  of  the  8192  x  8192  matrix  preserving  3-digit  accuracy  is  dune 
in  5  minutes  on  the  Sparcstation.  This  compares  to  95  days  (estimated) 
for  inverting  the  dense  matrix  by  Gauss-Jordan  and  to  24  minutes  for  one 
dense  matrix- vector  multiplication  of  that  size. 

The  condition  number  of  the  problem,  as  approximated  by  the  product  of  the 
row-sum  norms  of  V{1  —  T)UT  and  its  computed  inverse,  is  3  (independent  of 
size).  Five  iterations  were  required  by  the  Schulz  method  to  achieve  convergence. 
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Figure  5:  The  matrices  constructed  in  the  transformation  of  I  —  T.  matri¬ 
ces  Rq,...,R3  defined  in  Eq.  (30),  are  shown  for  kernel  I\(x,t)  =  log  j.r  —  t.\. 
e  —  10  ,  and  n  =  64.  The  matrix  R4  looks  like  R3  and  is  not  shown. 

In  Fig.  5  we  show  stages  in  the  transformation  of  the  matrix  I  —  T.  In 
particular,  for  t  =  10-3  and  n  =  64,  the  matrices  /?0,  . . . ,  Ri-i  defined  in  Eq.  (30) 
are  shown.  In  addition,  for  n  =  128  the  transformed  matrix  U(I  —  T)UT  and  its 
inverse  are  shown  in  Fig.  6. 

In  the  next  set  of  examples,  for  which  results  are  displayed  in  Table  2.  we 
used  the  wavelet  basis  of  order  k  =  8.  We  observe: 

1.  The  bandwidths  of  the  operator  matrix  and  its  inverse  are  less  for  k  =  8 
than  for  k  =  4.  The  inversion  times  are  correspondingly  smaller. 

2.  The  time  required  to  compute  the  operator  matrix  is  almost  four  times 
as  large  as  that  for  k  =  4.  This  is  due  to  the  cost  of  transforming  the 
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Figure  6:  Transformed  matrix  U(I  —  T)UT  (top)  and  its  inverse  (bottom)  arc 
shown  for  kernel  I\{x,  t )  —  log  jx  —  t\,  e  =  10-3,  and  n  =  128. 

near-diagonal  band,  which  is  twice  as  wide  for  k  =  8  as  for  k  =  4. 

3.  The  obtained  accuracy  exceeds  the  specified  precision  consistently. 
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Table  2:  The  integral  operator  1C  defined  by  the  formula  (/Cf)(x)  =  f(.r)~ 
fo  log  |x  —  t\f(t)dt  is  discretized,  transformed  to  the  wavelet  coordinates  with 
k  =  8,  and  inverted.  (See  Table  1  and  text.) 

Transform.  Inversion  Cfi 


£ 

n 

Ni 

t\ 

N  2 

*2 

Error 

IQ-2 

64 

5.8 

4 

6.2 

1 

0.161E-02 

128 

5.0 

16 

5.5 

2 

0.368E-02 

256 

3.3 

22 

3.6 

3 

0.184E-02 

512 

2.7 

46 

2.6 

4 

0.11.3E-02 

1624 

1.8 

62 

1.8 

4 

0.177E-02 

2648 

1.4 

182 

1.4 

5 

0.170E-02 

4666 

1.2 

363 

1.." 

8 

0.628E-03 

8162 

1.1 

726 

1.1 

11 

0.166E-02 

in-3 

64 

13.4 

5 

14.5 

8 

0.373E-03 

128 

14.2 

13 

15.5 

21 

0.332E-03 

256 

13.5 

28 

14.5 

46 

0.256E-03 

51" 

12.7 

57 

13.6 

90 

0.225E-03 

1624 

10.2 

114 

11.1 

134 

0.168E-03 

2048 

7.7 

221 

8.3 

176 

0.176E-03 

4666 

4.6 

426 

5.2 

185 

0.1 74E-03 

8162 

3.5 

818 

3.7 

208 

0.173E-03 

in-4 

64 

21.8 

6 

23.0 

23 

0.280E-04 

128 

26.3 

15 

28.0 

81 

0.253E-04 

256 

28.7 

35 

31.0 

235 

0.246E-04 

512 

28.4 

75 

30.6 

538 

0.184E-04 

1624 

25.5 

146 

27.2 

966 

0.625E-05 

2048 

22.0 

267 

23.8 

1736 

0.866E-05 

4666 

17.7 

561 

16.1 

2610 

0.768E-05 

4.  As  for  k  =  4 

,  the 

scaling  wi 

th  size 

n  is  li 

near  lor  the  transformation  step 

and  sublinear  for  the  inversion  step. 

In  the  final  set  of  examples  in  which  uncorrected  quadratures  were  used,  we 
perform  computations  for  k  =  4  and  e  =  10-3,  with  various  operator  kernels. 
Table  3  presents  the  results.  The  first  three  kernels  contain  singularities  of  the 
types  s(x)  =  log(x)  and  s(x)  =  xa  for  a  =  ±j,  and  are  nonsymmetric  and  noil- 
convolutional.  It  is  readily  seen  that  the  bandwidth  is  strongly  dependent  on  the 
type  of  singularity,  with  the  singularity  x“’^‘  producing  the  greatest  bandwidth. 
We  mention  also  that  this  particular  integral  equation  is  poorly-conditioned: 


Table  3:  The  integral  operator  A Z  defined  by  the  formula  ( Kf){x )  =  /(. r)  — 
/q1  I\(x,t)  f(t)  dt,  for  nonsymmetric,  nonconvolutional  kernels  K(x,t)  shown  be¬ 
low,  is  discretized,  transformed  to  the  wavelet  coordinates  with  k  =  4  and 
e  =  1 0~ 3 ,  and  inverted.  (See  Table  1  and  text.) 

Transform.  Inversion  C2 

I\(x,t)  n  N\  t\  N2  t2  Error 


cos(x<2)  log  |x  -  £| 

64 

18.2 

2 

20.2 

15 

0.31SE-03 

128 

18.6 

5 

20.4 

37 

0.302E-03 

256 

17.9 

11 

19.8 

82 

0.30  IE- 03 

512 

14.9 

22 

16.3 

131 

0.284E-03 

1024 

12.9 

42 

14.7 

242 

0.315E-03 

2048 

8.5 

76 

9.5 

283 

0.241E-03 

4096 

5.5 

137 

6.1 

291 

0.23  IE- 03 

8192 

3.6 

252 

4.3 

310 

0.230E-03 

cos(x<2)|x  -  t\~1/2 

64 

27.2 

3 

28.9 

32 

0.256E— 03 

128 

31.6 

7 

34.1 

122 

0.357E-03 

256 

35.6 

16 

40.6 

454 

0.434E-03 

512 

37.3 

35 

46.3 

1509 

0.643E-03 

1024 

34.5 

72 

45.4 

4166 

0.821E-03 

COS(xt2)|x  -  d1/2 

64 

6.8 

2 

7.3 

2 

0.303E-03 

128 

4.4 

4 

4.7 

2 

0.204E-03 

256 

2.9 

8 

3.0 

3 

0.209E-03 

512 

2.1 

15 

2.3 

3 

0.165E-03 

1024 

1.5 

30 

1.5 

3 

0.208E-03 

2048 

1.4 

60 

1.4 

6 

0.909E— 03 

4096 

1.1 

119 

1.2 

7 

0.614E-03 

8192 

1.1 

242 

1.1 

12 

0.666E-03 

(1  +  ^sin(100x))x 

64 

30.5 

3 

33.8 

44 

0.344E-03 

log|x  -  t\ 

128 

31.8 

6 

35.1 

103 

0.363E-03 

256 

21.2 

12 

24.1 

119 

0.34SE-03 

512 

18.6 

23 

20.7 

225 

0.372E-03 

1024 

15.8 

45 

18.4 

404 

0.392E— 03 

2048 

10.6 

82 

12.2 

466 

0.355E— 03 

4096 

6.4 

145 

7.4 

497 

0.336E— 03 

8192 

4.0 

265 

4.6 

510 

0.331E-03 

the  condition  numbers  of  the  discretizations  for  n  =  63, 128. 250. 512.  102  1  are 
9. 17, 34, 9S, 469,  respectively. 
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The  fourth  kernel  provides  an  example  with  an  oscillatory  coefficient  p(x)  = 
(1  +  |sin(100x)).  The  bases  developed  in  §3.4,  which  depend  on  p,  are  used  to 
transform  the  discretized  integral  operator  to  sparse  form.  We  see  in  Table  3 
that  the  inverse  is  also  very  sparse. 

5.2  Solution  of  Integral  Equations 

In  the  preceding  subsection,  we  examined  the  characteristics  of  various  integral 
operators  and  their  inverses  in  wavelet  coordinates.  We  used  completely  straight¬ 
forward  discretizations;  the  quadratures  represented  sums  of  the  integrands  at 
equispaced  points  (excluding  singular  points).  Such  simple  quadratures  converge 
too  slowly  to  the  integral  operators  to  be  of  much  use  in  solving  integral  equa¬ 
tions,  and  we  now  turn  to  the  high-order  quadratures  developed  in  [3]. 

We  first  present  examples  which  correspond  to  the  various  kernels  already 
tested  and  shown  in  Table  3.  In  Table  4  we  tabulate  the  results,  and  bandwidth 
differences  from  Table  3  reflect  the  effect  of  the  quadratures. 

For  the  remaining  examples  we  choose  integral  equations  that  can  be  solved 
analytically,  so  that  the  accuracy  of  the  method  can  be  checked.  We  consider  a 
class  of  integral  equations  with  logarithmic  kernel, 

f{x)-p{x)  [  log  |x  —  f|  f(t)  dt  =  gm(x),  x  €  [0, 1],  (32) 

Jo 

where  the  right  hand  side  gm  is  chosen  so  that  the  solution  /  is  given  by  the 
formula  /(x)  =  sin(mx).  The  integration  can  be  performed  explicitly,  yielding 

/  log  |x  —  t|  sin(?7it)  dt  =  log(x)  —  cos(m)  log(l  —  x) 

Jo 

—  cos(mx)[Ci(mx)  —  Ci(m(l  —  x))] 

—  sin(mx)[Si(mx)  +  Si(m(l  —  x ))] , 

where  Ci  and  Si  are  the  cosine  integral  and  sine  integral  (see,  e.g.,  [1],  p.  231). 
Equation  (32)  clearly  requires  quadratures  with  increasing  resolution  as  m  in¬ 
creases;  for  our  examples  we  let  n  =  m,  which  corresponds  to  2~  points  per 
oscillation  of  the  right  hand  side  gm. 

Initially  we  choose  coefficient  p(x)  =  1.  The  results  are  given  in  Table  5.  Here 
the  error  shown  is  the  error  of  the  computed  solution  relative  to  the  true  solution 
of  the  integral  equation.  Many  of  the  observations  of  the  preceding  examples  can 
be  repeated  here;  additionally,  we  make  the  following  comments: 

1.  The  bandwidths  are  greater  than  lor  the  uncorreued  quadratures,  but  this 
effect  generally  decreases  with  increasing  size. 
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Table  4:  The  integral  operator  fC  defined  by  the  formula  (ICf)(x)  =  /(.r)  — 
fo  K(x,t)  f(t)  dt,  for  nonsymmetric,  nonconvolutional  kernels  K(x,t)  shown  be¬ 
low ,  is  discretized  with  the  corrected  trapezoidal  rules,  transformed  to  the  wavelet 
coordinates  with  k  =  4  and  e  =  10-3,  and  inverted.  (Compare  to  Table  3.) 


Transform.  Inversion  C? 


N  i  ti _ N?  t2 _ Error 

28.3  4  31.6  38  0.164E-03 

31.5  9  34.3  103  0.162E-03 

30.8  21  33.9  221  0.172E-03 

27.0  41  29.7  370  0.177E-03 

21.0  80  23.7  454  0.357E-03 

14.8  143  17.2  566  0.317E-03 

9.5  250  10.4  555  0.282E-03 

5.8  448  6.9  665  0.271E-03 

cos(x<2)|x  -  t\~ll2  64  32.4  4  39.8  87  0.133E-02 

128  38.3  10  45.7  251  0.412E-03 

256  42.7  23  49.3  638  0.464E-03 

512  45.1  51  51.3  1494  0.562E-03 

1024  46.2  110  52.1  3309  0.635E-03 

cos(x<2)|x  -  <|1/2  64  10.4  3  18.4  9  0.867E-03 

128  7.6  6  13.8  13  0.526E-03 

256  5.1  13  9.3  16  0.358E-03 

512  3.3  25  5.2  15  0.292E-03 

1024  2.3  48  3.1  15  0.201E-03 

2048  1.9  96  2.3  20  0.393E-03 

4096  1.5  188  1.7  25  0.405E-03 

8192  1.3  374  1.4  36  0.404E-03 


2.  The  integral  equations  are  solved  to  within  the  specified  precision  in  every 
case  but  one.  The  exception,  for  e  =  10-4  and  n  =  64,  is  likely  due  to  the 
small  number  of  quadrature  points  and  high  specified  precision. 

3.  An  integral  equation  requiring  an  8192-point  discretization  is  solved  to  3- 
digit  accuracy  in  less  than  20  minutes  on  the  Sparcstation. 

For  our  second  set  of  integral  equations,  we  let  the  coefficient  p  be  the  os¬ 
cillatory  function  given  by  the  formula  p{x)  =  1  -f  |sin(100x).  We  carry  out 
the  transformation  described  in  §3.4  to  solve  the  integral  equation  32.  The  re¬ 
sults  are  shown  in  Table  6  and  as  with  Table  5  the  error  refers  to  the  error  of 
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Table  5:  The  integral  equations  f(x)  —  log  |i  —  t\  f(t)  dt  =  gm{x),  for  which  an 
explicit  solution  is  known,  are  solved  by  the  methods  of  this  chapter  (compare  to 
Table  1  and  see  text).  For  t  =  10— 2 , 10-3, 10-4  we  set  k  =  4,4,8,  respectively. 


=T 


Transform.  Inversion 


C? 


n,  m 

N  i 

ti 

n2 

<2 

Error 

64 

11.4 

3 

14.4 

1 

0.283E-02 

128 

10.7 

7 

13.2 

14 

0.212E-02 

266 

8.6 

13 

10.6 

20 

0.140E-02 

612 

6.4 

26 

7.6 

26 

0.112E-02 

1024 

3.6 

48 

4.6 

28 

0.821E-03 

2048 

1.0 

00 

2.3 

21 

0.032E-03 

4006 

1.3 

174 

1.6 

16 

0.674E— 03 

8102 

1.1 

344 

1.1 

13 

0.400E-03 

64 

27.7 

4 

31.3 

36 

0.236E-03 

128 

31.0 

0 

34.2 

00 

0.160E-03 

266 

30.6 

20 

33.6 

216 

0.161E— 03 

612 

27.6 

41 

30.2 

377 

0.130E— 03 

1024 

21.7 

70 

24.4 

470 

0.607E— 03 

2048 

16.6 

143 

18.1 

604 

0.470E— 03 

4006 

0.7 

248 

10.6 

670 

0.416E— 03 

8102 

6.0 

444 

7.3 

600 

0.364E— 03 

64 

37.2 

8 

46.0 

78 

0.127E-03 

128 

47.1 

23 

66.6 

278 

0.473E— 04 

266 

62.0 

64 

60.0 

746 

0.311E-04 

612 

66.0 

118 

61.4 

1701 

0.100E-04 

1024 

62.3 

248 

67.2 

3287 

0.734E-06 

10 


10 


-3 


10 


-4 


the  computed  solution  relative  to  the  true  solution  of  the  integral  equation.  For 
the  oscillatory  coefficient  we  see  performance  similar  to  the  constant-coefficient 
problem,  but  the  cost  is  higher. 

6  Generalizations  and  Applications 

In  this  paper,  we  have  constructed  a  new  class  of  vector-space  wavelet  bases  in 
which  a  variety  of  integral  operators  are  represented  as  sparse  matrices.  The 
inverses  of  these  matrices  are  also  sparse,  a  fact  which  enables  the  corresponding 
integral  equations  to  be  solved  rapidly.  We  have  asserted  that  the  time  complexity 
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Table  6:  The  integral  equations  f(x)—p(x)  log  |x  — 1\  f(t )  dt  =  gm(x),  for  which 
an  explicit  solution  is  known,  are  solved  by  the  methods  of  this  chapter  (compart 
to  Table  1  and  Sec  text).  For  t  ~  10-2, 10-3, 10-4  we  set  k  =  4,4,8,  respectively. 

Transform.  Inversion  C? 

e  n,m  N\  ti  N 2  t2  Error 

in-2 

64 

19.7 

4 

23.9 

18 

6.360E— 02 

128 

17.7 

8 

21.6 

36 

0.182E-02 

256 

12.6 

15  ■ 

14.6 

47 

0.174E— 62 

512 

8.4 

29 

9.8 

57 

0.112E-02 

• 

1624 

4.7 

55 

5.7 

56 

0.104E— 02 

2648 

2.4 

163 

2.7 

45 

0.902E-03 

4696 

1.6 

198 

1.7 

38 

0.720E-03 

8192 

1.3 

392 

1.3 

35 

0.543E-03 

10-3 

64 

36.2 

4 

41.3 

63 

6.22SE-02 

128 

46.8 

16 

47.6 

186 

0.209E-03 

256 

46.5 

23 

47.3 

427 

0.177E— 63 

512 

34.7 

46 

46.9 

712 

0.125E-03 

1924 

26.6 

87 

32.5 

1642 

0.134E— 03 

2648 

18.7 

158 

22.5 

1065 

0.5.97E— 03 

4696 

12.2 

281 

14.2 

1127 

0.529E-03 

8192 

7.2 

562 

8.4 

1164 

0.461E— 03 

in-4 

64 

47.6 

9 

58.2 

123 

0.236E— 62 

128 

66.7 

25 

77.3 

479 

0.186E-03 

256 

64.1 

59 

81.2 

1264 

0.124E— 03 

512 

62.5 

128 

76.3 

2492 

0.125E-04 

1624 

58.8 

267 

69.3 

4672 

0.862E-05 

for  an  n-point  discretization  is  bounded  by  order  0(n  log2  n),  but  observed  order 

0(n )  performance 

in  practice.  This  cost  should  be  contrasted  with  a  cost  of 

order  0(n 2)  for  direct  application  of  a  dense  matrix,  and  order  0(n3)  for  direct 

inversion. 

A  number  of  limitations  exist 

in  the  procedures  described  above.  These  re- 

strictions  may  be  c 

ategorized  as  “ 

software  limitations”  , 

and  “research  questions". 

We  discuss  software  limitations  first. 
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6.1  Software  Limitations 

Throughout  the  paper,  we  have  assumed  that  the  size  of  the  problem  n  has  the 
form  n  =  2lk  for  some  /.  This  restriction  is  not  fundamental;  it  merely  simplifies 
the  software. 

A  second  software  restriction  is  the  assumption  of  only  diagonal  singularities. 
This  case  is  an  important  one  in  practice,  but  in  certain  situations  one  may 
encounter  singularities  or  near-singularities  off  the  main  diagonal.  The  scheme 
described  in  §4.2  for  transformation  of  a  matrix  to  wavelet  bases  can  be  readily 
revised  to  an  adaptive  scheme,  which  works  as  follows:  an  m  x  m  submatrix 
A  is  transformed  to  wavelet  coordinates,  under  the  assumption  that  it  can  be 
approximated  to  high  precision  in  each  direction  by  a  polynomial  of  degree  less 
than  k.  This  assumption  is  then  checked  by  dividing  A  into  four  submatrices,  each 
of  dimension  m/2  x  m/2,  transforming  each  submatrix,  and  “glueing”  the  pieces 
together.  If  the  results  from  the  two  computations  match  (to  high  precision),  no 
further  refinement  of  the  original  submatrix  is  needed.  Otherwise,  the  procedure 
is  repeated  recursively  on  the  m/2  x  m/2  submatrices.  The  cost  of  this  adaptive 
procedure  is  roughly  5  times  as  great  as  the  cost  of  a  static  procedure  in  which 
the  structure  of  the  singularities  is  known  a  priori. 

6.2  Research  Questions 

The  list  of  research  issues  is  of  course  much  longer.  One  of  the  most  pressing  issues 
is  the  generalization  to  two  and  three  dimensions.  Although  conceptually  the 
generalization  of  the  wavelet  bases  to  several  dimensions  is  quite  straightforward 
(see,  e.g.,  [2]),  actual  procedures  to  perform  the  required  orthogonalizations  have 
not  been  developed.  Also,  the  issue  of  high-order  quadratures  for  two  and  three 
dimensions  has  not  been  resolved. 

Another  question  is  whether  similar  “custom-constructed”  bases  can  be  used 
to  create  sparse  representations  of  integral  operators  with  oscillatory  kernels. 
Initial  efforts  in  this  direction  for  a  limited  class  of  such  operators,  in  particular  for 
Fourier  transforms  with  non-equispaced  points  and  frequencies,  appear  promising 

[9]- 

6.3  Applications 

In  this  paper  the  primary  application  of  our  new  wavelet  bases  has  been  the 
solution  of  second-kind  integral  equations.  The  bases  are  very  effective  for  the 
fast  solution  of  a  wide  class  of  such  problems.  In  addition,  many  other  classes  of 
problems  can  be  solved  efficiently  using  these  techniques.  We  list  a  few  of  th^se 
problem  types. 
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1.  Elliptic  partial  differential  equations  rewritten  as  integral  equations  by 
the  Lippman-Schwinger  method,  in  which  the  Green’s  functions  are  non- 
oscillatory. 

2.  Evolution  of  homogeneous  parabolic  PDEs  with  constant  or  periodic  bound¬ 
ary  conditions,  by  explicit  time  steps.  This  method  consists  of  repeated 
squarings  of  the  operator  for  a  single  time  step,  leading  to  an  order  0{n  log  t ) 
algorithm  for  evolving  an  n- point  discretization  for  t  time  steps. 

3.  Evolution  of  general  parabolic  PDEs  by  implicit  time  steps,  in  which  the 
elliptic  problem  on  each  time  step  is  solved  in  wavelet  coordinates. 

4.  Evolution  of  hyperbolic  PDEs  by  a  method  of  operator  squaring  analogous 
to  the  scheme  proposed  for  homogeneous  parabolic  PDEs  above. 

5.  Problems  of  potential  theory  and  pseudo-differential  operators. 

6.  Signal  compression,  including  signals  of  seismic,  visual,  and  vocal  origin. 
There  is  also  reason  to  expect  that  analysis  of  such  compressed  data  will 
be  simpler  than  analysis  of  data  resulting  from  less  efficient  compression 
schemes. 

In  this  paper  we  strayed  from  the  original  mathematical  definition  of  wavelets 
to  construct  classes  of  bases  tailored  for  numerical  computation.  The  basis  vec¬ 
tors’  principal  properties  of  local  support  and  vanishing  moments  lead  to  sparse 
representations  of  functions  and  operators  that  are  smooth  except  at  a  small 
number  of  singularities.  There  is  little  doubt  that  other  bases  can  be  constructed 
along  similar  lines  to  possess  various  properties.  One  current  challenge  is  the 
construction  of  bases  suitable  for  the  efficient  representation  of  a  variety  of  oscil¬ 
latory  operators. 
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